#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <pthread.h>
#define NITERATION 23

/*Struct that contains the values of An, Bn and Pn*/
typedef struct{
    mpf_t a;
    mpf_t b;
    mpf_t p;
}Data;

/*Global Variables*/
Data D[NITERATION];
pthread_mutex_t m1[NITERATION], m2[NITERATION], m3[NITERATION];

/*THe inicial functions for each thread*/
void *Value_A(void *k){
    int i;
    mpf_t v_sqrt;
    mpf_t temp;

    mpf_init_set_ui(v_sqrt, 0l);
    mpf_init_set_ui(temp, 0L);

    for(i=0; i<NITERATION-1; i++){
        //Update of An+1
        mpf_sqrt(v_sqrt, D[i].a);

        mpf_ui_div(temp, 1L, v_sqrt); // temp = 1 / sqrt(A(n))
        mpf_add(D[i+1].a, v_sqrt, temp); // A(n+1) = sqrt(A(n)) + (1/sqrt( A(n)))
        mpf_div_ui(D[i+1].a, D[i+1].a, 2L); // A(n+1) = (sqrt(A(n)) + (1/sqrt(A(n)))) / 2

        pthread_mutex_unlock(&m1[i]);
        pthread_mutex_unlock(&m2[i]);
    }
    pthread_exit(NULL);
}

void *Value_B(void *k){
    int i;
    mpf_t v_sqrt;
    mpf_t temp;
    mpf_init_set_ui(v_sqrt, 0l);
    mpf_init_set_ui(temp, 0L);

    for(i=0; i<NITERATION-1; i++){

        pthread_mutex_lock(&m1[i]);

        mpf_sqrt(v_sqrt, D[i].a);
        //Update of Bn+1
        mpf_add_ui(temp, D[i].b, 1L); //temp = 1 + b(n)
        mpf_mul(D[i+1].b, temp, v_sqrt); //b(n+1) = (1 + b(n)) * sqrt(A(n))
        mpf_add(temp, D[i].a, D[i].b); //temp = a(n) + b(n)
        mpf_div(D[i+1].b, D[i+1].b, temp); // B(n+1) = (1 + b(n)) * sqrt(A(n)) / a(n) + b(n)
        pthread_mutex_unlock(&m3[i]);
    }
    pthread_exit(NULL);
}

void *Value_P(void *k){

    int i;
    mpf_t temp;

    mpf_init_set_ui(temp, 0L);
    sleep(1);

    for(i=0; i<NITERATION-1; i++){
        pthread_mutex_lock(&m2[i]);
        pthread_mutex_lock(&m3[i]);

        //Update of Pn+1
        mpf_add_ui(temp, D[i+1].a, 1L); //temp = A(n+1) + 1
        mpf_mul(temp, temp, D[i].p); // temp = (A(n+1) + 1) * p(n)
        mpf_mul(D[i+1].p, temp, D[i+1].b); // P(n+1) = ( (A(n+1) + 1) * p(n) ) * b(n+1)
        mpf_add_ui(temp, D[i+1].b, 1L); // temp = (1 + B(n+1) )
        mpf_div(D[i+1].p, D[i+1].p, temp); // P(n+1) = ( (A(n+1) + 1) * p(n) ) * b(n+1) / (1 + B(n+1) )
    }
    pthread_exit(NULL);
}

int main(){
    /*Declaration and setting the inicial values on Variables*/
    pthread_t Thread_a, Thread_b, Thread_p;
    int i;
    int flag;
    mpf_t v_sqrt;
    mpf_t temp;

    for(i=0; i<NITERATION; i++) pthread_mutex_init(&m1[i], NULL);
    for(i=0; i<NITERATION; i++) pthread_mutex_lock(&m1[i]);

    for(i=0; i<NITERATION; i++) pthread_mutex_init(&m2[i], NULL);
    for(i=0; i<NITERATION; i++) pthread_mutex_lock(&m2[i]);

    for(i=0; i<NITERATION; i++) pthread_mutex_init(&m3[i], NULL);
    for(i=0; i<NITERATION; i++) pthread_mutex_lock(&m3[i]);

    mpf_set_default_prec(40000016);

    for(i=0; i<NITERATION; i++){
        mpf_init_set_ui(D[i].a, 0L);
        mpf_init_set_ui(D[i].b, 0L);
        mpf_init_set_ui(D[i].p, 0L);
    }

    mpf_sqrt_ui(D[0].a, 2L);

    mpf_add_ui(D[0].p, D[0].a, 2L);

    /*Creation of the Threads. Each one execute a different function */
    flag = pthread_create(&Thread_a, NULL, Value_A, NULL);
    if(flag) exit(-1);

    flag = pthread_create(&Thread_b, NULL, Value_B, NULL);
    if(flag) exit(-1);

    flag = pthread_create(&Thread_p, NULL, Value_P, NULL);
    if(flag) exit(-1);

    /*Waits for all the threads to continue the program*/
    pthread_join(Thread_a, NULL);
    pthread_join(Thread_b, NULL);
    pthread_join(Thread_p, NULL);

    /*The number of pi is printed in the standart output*/
    mpf_out_str(stdout, 10, 40000016, D[NITERATION-1].p);

    /*Clear the variables*/
    for(i=0; i<NITERATION-1; i++){
        mpf_clear(D[i].a);
        mpf_clear(D[i].b);
        mpf_clear(D[i].p);
    }

    return 0;
}
